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Objectives  and  Approach 

The  seismological  problem  most  relevant  to  nuclear  monitoring  under  a  Comprehensive  Test 
Ban  Treaty  (CTBT)  is  to  discriminate  small  nuclear  explosions  from  chemical  explosions  and 
small  earthquakes  in  a  variety  of  geological  settings.  The  dual  challenges  posed  by  the  desire  for 
a  small  detection  threshold  and  a  global  reach,  which  includes  areas  of  low  seismicity,  call  for  an 
approach  to  discriminate  nuclear  explosions  from  chemical  explosions  and  earthquakes  that 
simultaneously  uses  multiple  diagnostics,  such  as  Ms-mb  ratio,  radiation  pattern,  and  amplitude 
spectrum.  A  general  method  for  studying  seismic  sources  that  encompasses  many  traditional 
source  parameter  discriminants  is  based  on  the  recovery  of  the  broad-band  moment-rate  tensor, 
M(m).  This  moment  tensor  can  be  modeled  as  the  product  of  a  complex  function,  /(G)) ,  whose 
inverse  Fourier  transform  is  the  source-time  function,  f(t),  and  a  time-independent  mechanism 
tensor,  M,  usually  taken  to  be  the  (zeroth-order)  spatial  moment  of  the  stress  glut-rate  tensor 
[Backus,  1977;  Silver  and  Jordan,  1982].  In  explosion  seismology,  the  source-time  function  is 
often  equated  to  the  reduced  displacement  potential  (RDP),  y{t),  via  f(t)  =  Aitpa2  y/(t)  [e.g., 
Denny  and  Johnson,  1991].  More  generally,  however,  several  processes  contribute  to  M(g>), 
these  include  the  release  of  tectonic  prestress,  spall  and  slapdown,  as  well  as  the  RDP  of  the 
explosion  itself. 

The  comprehensive  nature  of  an  approach  based  on  recovering  M(g>)  is  appealing.  For 
example,  writing  /(G))  =  A{(0)em<0)  in  polar  form  illustrates  that  the  amplitude  spectrum, 
A{(o),  is  just  the  modulus  of  /(G)).  Furthermore,  to  the  extent  that  good  synthetic  seismograms 
can  be  computed  from  M(o) ,  relative  excitations  of  different  waves  (e.g.,  Ms-mb  ratio)  can  also 
be  obtained.  In  addition,  several  new  types  of  discriminants  are  potentially  derivable  from 
/(G)) .  The  moment-tensor  approach  is  amenable  for  setting  up  inverse  problems  to  separate  the 
isotropic  explosion  component,  M/(g))  =  M/(g>)I,  from  deviatoric  sources  of  wave  excitation, 
Md(g)).  The  shallow  depth  of  explosion  sources  makes  an  inversion  for  all  six  moment-tensor 
components  (M/(G))  and  the  five  independent  elements  of  M £>(&)))  difficult.  Long-period 
surface  waves,  for  example,  are  able  to  constrain  only  three  components  [Given  and  Mellman, 
1986;  Ekstrom  and  Richards,  1994].  In  addition  to  mr6  and  mr<t> ,  poorly  constrained  since  their 
Green's  functions  vanish  for  surface  sources,  an  independent  estimation  of  mrr  and  mee  + 
is  difficult  using  only  surface  waves  [Patton,  1988].  To  overcome  these  limitations  Patton 
[1988;  1991]  performed  moment-tensor  inversions  of  regional  observations  from  NTS 
explosions  using  both  fundamental-mode  and  higher-mode  surface  waves.  The  higher-mode  data 
provide  additional  constraints  on  the  source,  since  Green's  functions  of  mrr  and  meQ  + 
differ  substantially  even  for  shallow  sources  for  these  waves. 
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An  important  goal  of  this  study  is  to  utilize  our  broad-band  source  recovery  algorithms 
employing  both  surface- wave  and  body-wave  data  to  determine  M/(o))  and  M £,(&))  and  to 
quantify  the  different  contributions  to  M D{(0),  especially  the  tectonic  release  component, 
Mr(£0),  whose  contribution  dominates  the  deviatoric  component  at  low  frequencies.  Tectonic 
release  contributes  to  source  parameter  estimates,  such  as  the  surface-wave  magnitude  Ms  [e.g., 
Patton ,  1991],  and  has  the  potential  to  bias  both  discriminants  and  minimum-yield  detection 
maps,  in  particular  when  the  amount  of  tectonic  release  is  significant.  A  quantitative  measure  for 
the  tectonic  component  is  the  "F  factor",  defined  by  Toksoz  et  al.  [1964]  as  the  ratio  of  the 
Rayleigh-wave  amplitude  generated  by  the  explosion  to  that  generated  by  the  tectonic  release. 
For  a  Poisson  solid  Mt/Mj  =2/3 F  [Wallace,  1991]. 

The  presence  of  a  tectonic  release  component  in  seismic  recordings  from  nuclear  explosions, 
manifested  in  the  presence  of  S/Hype  seismic  waves,  has  been  recognized  since  the  1960s  [Press 
and  Archambeau,  1962;  Brune  and  Pomeroy,  1963].  Models  proposed  for  explaining  the 
generation  of  Sh- type  energy  include  explosion-triggered  slip  at  preexisting,  nearby  faults  [Aki 
and  Tsai,  1972],  release  of  tectonic  prestress  around  the  explosion  source  [Archambeau,  1972], 
and  a  combination  of  the  two  [Wallace,  1991].  Important  unresolved  questions  associated  with 
the  tectonic  release  component  are:  (1)  What  is  the  scaling  of  the  tectonic  component  with 
frequency?  (2)  Is  there  a  change  in  the  character  of  tectonic  release  with  explosion  yield/depth? 
(3)  How  does  tectonic  release  vary  in  different  geologic  settings? 

Most  previous  studies  of  tectonic  release  have  focused  on  the  Nevada  Test  Site  [Aki  and  Tsai, 
1972;  Wallace  et  al.,  1983;  Wallace  et  al.,  1985;  Given  and  Mellman,  1986;  Patton,  1988]  and 
the  Shagan  River,  Kazakhstan  nuclear  test  site  [Helle  and  Rygg,  1984;  Given  and  Mellman, 
1986;  Walter  and  Patton,  1990;  Ekstrom  and  Richards,  1994],  while  little  attention  has  been 
devoted  to  investigating  the  tectonic  component  from  nuclear  tests  near  Lop  Nor,  China  [Zhang, 
1996]. 

Furthermore,  most  studies  have  utilized  either  surface  waves  [e.g..  Given  and  Mellman,  1986; 
Ekstrom  and  Richards,  1994;  Zhang,  1996]  or  body  waves  [Wallace  et  al.,  1983;  Wallace  et  al., 
1985]  alone,  making  a  broad-band  estimate  of  M(n>)  difficult.  Therefore,  the  variation  of  F  with 
frequency  has  received  little  attention,  although  scarcity  of  high-frequency  Sh  waves  indicates 
that  it  decreases  with  frequency.  Our  approach  based  on  a  broad-band  source  recovery  will  allow 
us  to  constrain  this  quantity  and  improve  our  physical  understanding  of  explosion  source 
processes. 

Another  important  question  concerns  the  scaling  of  tectonic  release  for  different  explosions. 
Some  authors  have  suggested,  for  example,  that  F  decreases  with  time  for  explosions  detonated 
in  close  proximity  to  each  other  [Aki  and  Tsai,  1972;  Wallace,  1991],  and  that  it  is  affected  by 
other  factors  such  as  explosion  yield  and  depth  of  burial  [Aki  and  Tsai,  1972].  From  a  moment- 
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tensor  inversion  of  16  explosions  on  Pahute  Mesa,  Patton  [1991]  concluded  that  the  character  of 
the  tectonic  component  changes  with  explosion  yield,  going  from  vertical  strike-slip  mechanisms 
to  dip-slip  reverse  faulting  at  yields  of  about  200-300  kT.  Since  Ms  is  biased  high  (low)  for 
strike-slip  (thrust)  mechanisms  [Patton,  1980],  it  is  important  to  test  this  hypothesis  at  other 
nuclear  test  sites. 

In  this  study,  we  consider  explosions  from  the  Chinese  nuclear  test  site  at  Lop  Nor.  We  have 
collected  a  large  number  of  broad-band  and  long-period  three-component  seismograms  at 
regional  and  teleseismic  distances  for  several  nuclear  tests  from  1990-1996  spanning  a  range  of 
sizes.  Applying  our  broad-band  source  recovery  algorithms  to  this  data  set,  we  will  be  able  to 
improve  our  physical  understanding  of  the  processes  affecting  seismic  wave  excitation  from 
explosion  sources.  Gaining  new  insights  into  how  tectonic  release  scales  with  frequency  and 
explosion  yield/depth  should  also  lead  to  a  better  understanding  on  how  high-frequency 
discriminants  such  as  P/Lg  are  affected  by  tectonic  release. 

Recovering  broad-band  estimates  of  the  source-time  function  from  seismic  data  at  regional 
and  teleseismic  distances  is  a  difficult  inverse  problem,  requiring  precise  corrections  for 
deterministic  propagation  effects  and  careful  averaging  procedures  to  reduce  stochastic 
fluctuations  due  to  non-deterministic  scattering.  Our  approach  is  based  on  techniques  developed 
by  Jordan  and  coworkers  for  determining  earthquake  source-time  functions.  In  particular,  we 
adapt  waveform  isolation  and  analysis  tools  developed  by  Gee  and  Jordan  [1992]  and  Ihmle  and 
Jordan  [1995],  which  use  synthetic  seismograms  to  improve  the  localization  of  signal 
measurements  in  both  time  and  frequency  domains.  We  proceed  as  follows: 

(1)  For  each  observed  seismogram  Sj(t)  (j  =  station  index),  we  generate  six  Green  s  function 
seismograms  gj(t)  (k  =  component  index;  k  e  [rr,  66,  0</>,  r6  ,  r<p,  6<j>])  from  the 
best  available  structural  model.  At  teleseismic  distances,  we  use  a  global  tomographic 
model,  such  as  the  Preliminary  Reference  Earth  Model  (PREM)  of  Dziewonski  and 
Anderson  [1981]  in  conjunction  with  the  degree-12  structure  S12_WM13  of  Su  et  al. 
[1994],  Using  Green's  functions  calculated  for  an  aspherical  earth  model  reduces 
contamination  due  to  propagation  effects  and  enables  us  to  obtain  a  clearer  image  of  the 
source.  For  the  work  at  regional  distances,  we  adopt  appropriate  path-averaged  models 
[e.g.,  Gaherty  and  Jordan,  1995],  where  available.  We  use  an  iterative  inversion 
algorithm  to  determine  a  source  mechanism,  m(ft)),  from  a  series  of  narrow-band  filtered 
cross-correlations  between  Sj(t)  and  g){t).  (m(G>)  is  the  isomorphic  equivalent  of  the 
symmetric  tensor  M(a>)  [Silver  and  Jordan,  1982].) 
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(2)  Using  an  average  source  mechanism,  m  =  m(ot>o)>  we  generate  synthetic  seismograms 
Sj  ( t )  =  g j  (r)  •  m  assuming  that  the  source  is  a  step-function  at  r0  • 

(3)  Using  the  isolation-filtering  procedures  detailed  by  Gee  and  Jordan  [1992],  we  extract 
two  functions  of  frequency,  an  amplitude  ratio  Aj  (<u)  /  Aj(co)  and  a  differential  phase 
delay  Atj(co)  =  tj(co)-tj(co),  from  individual  waveforms  on  Sj(t).  At  teleseismic 
distances,  the  waveforms  we  analyze  typically  include  the  P  and  S  body  phases  and  the 
R{  and  G!  surface  waves.  Corrections  are  applied  to  account  for  the  amplitude  and  phase 
distortions  due  to  the  time-  and  frequency-localization  operations  intrinsic  to  the  analysis 
[Gee  and  Jordan,  1992]. 

(4)  We  average  Atj(O))  and  the  logarithm  of  Aj(co)/  Aj(co)  over  the  network  from  which 
we  obtain  estimates  of  At(co)  and  A(a>).  Note  that  averaging  is  not  done  over  the 
complex  spectrum,  f(co );  if  it  were,  the  small  phase  delays  associated  with  unmodeled 
propagation  effects  would  cause  destructive  interference  among  the  signals,  and  therefore 
lead  to  amplitude  bias. 

(5)  We  invert  the  estimates  of  A(co)  and  At(co)  for  fit),  subject  to  an  appropriate  set  of 
smoothing  constraints  and  functional  bounds. 

Research  Accomplishments 

Modifications  to  the  Analysis  Procedures.  We  have  adapted  our  algorithms  to  synthesize  broad¬ 
band  amplitude  and  phase-delay  estimates  from  regional  and  teleseismic  data,  developed  in  the 
context  of  our  research  on  earthquake  source  mechanisms  [Ihmle  et  al.,  1993;  Ihmle  and  Jordan, 
1994],  to  inversions  of  explosion  sources.  In  our  previous  earthquake  studies,  we  have  assumed 
that  the  source-time  function  is  non-negative.  Physically,  this  is  a  statement  that  there  is  no 
backslip  in  the  fault.  The  mathematical  advantage  of  this  assumption  is  that  a  lower  bound  is 
imposed  on  the  inversion  for  fit),  which  constrains  the  solution  by  reducing  the  size  of  the 
model  space  [Ihmle  and  Jordan,  1994].  This  assumption  is  not  generally  valid  for  nuclear 
explosions,  because  many  of  them  are  observed  to  have  overshoot,  which  means  that  fit) 
becomes  negative  over  some  time  interval  of  the  source.  However,  when  inverting  for  the 
deviatoric  (tectonic)  component  of  the  source-time  function,  foit),  non- negativity  should  be  a 
good  first-order  approximation.  Another  constraint  commonly  applied  to  earthquake  moment 
tensor  inversions  is  to  require  the  moment  tensor  to  be  deviatoric.  Such  a  constraint  is  clearly 
inappropriate  when  inverting  for  a  moment  tensor  that  is  composed  of  both  explosion  and 
deviatoric  components.  A  further  difficulty  encountered  when  inverting  for  explosion  source 
mechanisms  is  due  to  the  fact  that  the  depth  of  moment  release  is  very  close  to  the  earth's 
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surface.  At  shallow  depths  the  Green's  functions  gr0(t)  and  gr^(t)  almost  vanish,  since  the 
normal  component  of  the  stress  tensor  is  zero  at  the  free  surface.  This  means  that  mrQ  and  mr<p 
cannot  be  well  constrained.  A  common  practice  in  source  inversions  of  nuclear  explosions  and 
shallow  earthquakes  is  to  simply  set  these  two  components  equal  to  zero  [Kanamori  and  Given, 
1981;  Given  and  Mellman,  1986].  However,  this  procedure  is  adequate  only  for  narrow -band 
moment-tensor  inversions,  since  the  sensitivity  to  shallow  structure  becomes  larger  with 
increasing  frequencies,  which  leads  to  high-frequency  moment-tensor  solutions  that  are  biased 
with  respect  to  the  low-frequency  inversions  when  mr6  =  mr(p  =  0  is  assumed.  We  overcome 
this  problem  by  recasting  the  inverse  problem  using  a  projection  operator,  P*  =  I  -  gs  gj .  This 
operator  annihilates  any  possible  information  on  mrQ  and  mr(p  contained  in  the  Green's  functions 
g  or  in  the  data  s.  This  modifies  the  inverse  problem  from: 

g-m  =  gd  ■  [mrr , m96 , ,md<p]  +  gs- [mr6 , mr4> ]  =  s  (1) 

to 

Psg  m  =  Pss, 

and  allows  a  recovery  of  the  moment  tensor  over  a  broad  frequency  range.  Another  advantage  of 
the  projection  operator  approach  is  that  it  allows  a  self-consistent  error  analysis,  since  the 
posterior  model  covariance  matrix  will  reflect  the  constraints  used  in  obtaining  the  solution. 

Our  formalism  for  obtaining  the  source  mechanism  is  based  on  determining  the  elements  of 
the  moment  tensor  directly.  This  contrasts  with  methodologies  that  formulate  the  Green's 
functions  using  three  fundamental  faults  [Kanamori  and  Given,  1981]  plus  an  explosion  source 
[Given  and  Mellman ,  1986].  We  also  do  not  impose  any  a  priori  constraint  on  the  fault  geometry 
or  even  require  the  deviatoric  component  to  be  a  double  couple. 

Synthetic  Test  Case.  To  examine  the  performance  of  our  algorithms  we  construct  the  following 
test  case.  Using  a  fictitious  source  located  at  the  Lop  Nor  test  site  with  a  moment  tensor 
composed  of  an  isotropic  part  plus  a  deviatoric  mechanism  with  identical  scalar  moment,  we 
calculate  3-D  synthetic  seismograms  for  a  PREM  plus  S12_WM13  earth  structure  for  20  IRIS 
and  Geoscope  stations  at  regional  and  teleseismic  distances  to  simulate  the  data.  The  aspherical 
synthetics  are  calculated  employing  the  asymptotic  formalism  of  Woodhou.se  and  Dziewonski 
[1984].  The  deviatoric  moment  tensor  is  the  Harvard-CMT  solution  of  the  94/9/7  event  in 
Southern  Xinjiang,  China  (a  predominantly  strike-slip  earthquake  with  a  small  CLVD 
component  and  a  scalar  moment  of  7.7  x  1016  Nm)  mimicking  the  tectonic  release  component. 
Green’s  functions  are  calculated  for  the  1-D  PREM  model.  Our  imprecise  knowledge  of  3-D 
earth  structure  is  a  primary  source  of  signal-generated  noise  which  makes  an  inversion  for  the 
source-time  function  difficult.  Using  3-D  synthetics  as  data  and  1-D  Green's  functions  provides 
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a  proxy  for  difficulties  encountered  with  real  data.  Since  the  explosion  source  is  located  at 
shallow  depth  (1  km),  the  mr6  and  mr<p  components  of  the  moment  tensor  are  ill-constrained;  we 
therefore  project  out  these  components  using  the  algorithm  outlined  above.  Utilizing  the  fact  that 
transverse-component  (LHT)  seismograms  do  not  contain  any  energy  from  a  purely  isotropic 
source  our  moment-tensor  inversion  proceeds  in  two  steps.  First,  we  invert  the  LHT 
seismograms  and  obtain  an  estimate  for  the  tectonic  release.  Second,  we  invert  the  vertical- 
component  data,  which  yields  a  moment  tensor  that  is  a  combination  of  the  isotropic  component 
plus  the  tectonic  release. 

Figure  1  shows  the  results  of  the  initial  moment-tensor  inversion  using  long-period,  minor- 
arc  Love  waves  recorded  on  the  transverse  component.  We  apply  a  series  of  narrow-band  filters 
to  the  broad-band  cross-correlation  function  and  obtain  independent  moment-tensor  solutions 
for  each  frequency  band.  The  source  mechanism  plots  are  constructed  using  the 
methodology  described  by  Riedesel  and  Jordan  [1989].  Column  A  displays  the  three  principal 
axes  of  the  source  mechanism  (T,  N,  P)  together  with  the  unit  mechanism  vector  (X)  and  the 
canonical  unit  vectors  d,  1,  and  i,  representing  a  pure  double  couple,  a  compensated  linear 
vector  dipole,  and  a  pure  dilatation,  respectively.  The  contours  of  the  P-wave  first  motion  are 
shown  in  column  B.  This  approach  provides  a  comprehensive  visualization  of  the  source 
mechanism.  In  particular,  it  is  possible  to  directly  assess  the  uncertainties  in  the  source 
mechanism  and  quantify  the  contribution  of  an  isotropic  component,  a  very  useful  feature  when 
investigating  non-deviatoric  sources  (see  below).  While  all  solutions  are  strike-slip  mechanisms, 
the  source  mechanisms  vary  rapidly  between  adjacent  frequency  bands  and  differ  from  the 
projected  ground-truth  mechanism  (bottom  right).  This  instability  in  the  inversion  is  primarily  a 
result  of  mapping  phase-shifts  due  to  unmodeled  lateral  heterogeneities  along  the  propagation 
path  into  the  source  mechanism.  To  account  for  these  propagation  effects  we  iterate  our 
inversion  by  (1)  calculating  the  phase-shifts  between  data  and  synthetics  (computed  from  a 
moment  tensor  obtained  in  the  initial  inversion),  (2)  phase-shifting  the  Green  s  functions,  and  (3) 
reinverting  for  MD(co).  The  results  depicted  in  Figure  2  document  the  success  of  this 
methodology.  The  agreement  between  different  frequency  bands  and  with  the  projected  ground- 
truth  moment  tensor  is  excellent. 

Horizontally  polarized  shear  waves  provide  additional  constraints  on  (&>)  and  are  less 
contaminated  by  shallow  heterogeneity  than  surface  waves.  Figure  3  shows  the  results  of  the 
moment-tensor  inversion  using  Sh  waves,  adjusted  for  the  phase-shifts  due  to  aspherical 
heterogeneities.  The  agreement  with  the  Love-wave  solutions  is  good,  and  indicates  that  the  two 
data  sets  can  be  used  to  cross-check  the  inversion  results  and  provide  a  broad  frequency  band. 
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Figure  1.  Source  mechanism  plots  (A)  and  initial  P- wave  radiation  patterns  (B)  for  MD(e>)  in  2 
mHz  frequency  bands  ranging  from  9  to  31  mHz  plotted  using  equal  area  projection.  On  each 
plot  (A),  the  three  principal  axes  elf  e2,  e3  (eigenvectors  of  MD(co))  corresponding  to 
eigenvalues  A,  >  A2  ^  A3,  are  plotted  as  T,  N,  and  P  on  the  lower  focal  hemisphere.  The 
mechanism  is  characterized  by  the  unit  vector  A  =  XAtej,  which  can  be  compared  with  the 
canonical  unit  vectors  d,  1,  and  i,  representing  a  pure  double  couple,  a  pure  compensated  linear 
vector  dipole,  and  a  pure  dilatation,  respectively.  The  dashed  line  shows  the  locus  of  all 
deviatoric  mechanisms.  Vectors  are  plotted  as  downward  pointing  filled  triangles  if  they  are  on 
the  lower  focal  hemisphere  and  as  upward  pointing  if  they  are  projected  from  the  upper 
hemisphere.  For  a  pure  double  couple  the  axes  T,  N,  and  P  are  the  tension,  neutral,  and 
compression  axes,  respectively.  The  contours  of  the  P- wave  first  motion  are  shown  in  (B).  Solid 
(dashed)  contours  denote  compressional  (dilatational)  first  motion,  thick  contours  are  nodal  lines. 
Initial  moment-tensor  inversion  results  using  first-orbit,  transverse-component  Love  wave-forms 
from  a  synthetic  test  case.  "Data"  seismograms  are  calculated  using  PREM  plus  the  3D  model 
S12_WM13.  Green’s  functions  are  calculated  for  the  ID  PREM  model.  The  mrB  and  mr(p 
components  are  projected  out  from  both  solution  and  data  space.  The  moment  tensor  is 
constrained  to  be  deviatoric.  For  comparison,  the  projected  input  moment  tensor  ("Ground 
Truth")  is  shown  on  the  lower  right. 
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Figure  2.  Mechanism  plots  (A)  and  P-wave  radiation  patterns  (B)  from  the  Love-wave  inversion 
(phase- aligned)  for  the  test  case  described  in  Figure  1.  Prior  to  performing  the  moment-tensor 
inversions,  a  phase-shift,  A tj(co),  is  calculated  for  each  seismogram  in  each  frequency  band 
using  a  preliminary  Mc  from  the  initial  moment-tensor  inversion  (Figure  1).  The  Green's 
functions  are  phase-shifted  by  A tj(co)  and  a  new  inversion  for  M £>(&))  is  performed.  The 
agreement  between  the  frequency  bands  and  with  the  projected  ground-truth  moment  tensor  is 
remarkable. 
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Figure  3.  Mechanism  plots 
(A)  and  P-wave  radiation 
patterns  (B)  from  the  mo¬ 
ment-tensor  inversion  (after 
phase  alignment)  using  S# 
waveforms  for  the  same  test 
case  described  in  Figure  1. 
The  mrQ  and  mr(p  compo¬ 
nents  are  projected  out  and 
the  moment  tensor  is  con¬ 
strained  to  be  deviatoric. 
Results  from  the  body-wave 
inversions  agree  very  well 
with  the  surface- wave  results 
displayed  in  Figure  2. 


Network-averaged  amplitudes,  A{(o)l  A(co),  and  time-shifts,  Ar(cw),  for  the  Sh  waves 
obtained  from  cross-correlations  of  data  and  synthetics  (using  an  average  moment  tensor  from 
Figure  3)  are  displayed  in  Figure  4.  The  small  residuals  ( A(co)  /  A(co)  -  1  and  At(co)  ~  0)  are 
further  proof  of  the  success  of  the  moment-tensor  inversion.  The  residual  time  shifts  in  Figure  4, 
are  due  to  3D  propagation  effects  not  removed  by  the  network  average.  Using  these  amplitude 
and  time  shift  spectra  we  can  invert  for  the  source-time  function. 


Figure  4.  Network-averaged  amplitude  ratios,  A(co)/A(co),  and  time-shifts,  A t(co),  measured 
from  narrow-band  cross-correlations  of  Sh  waveforms  using  an  average  moment  tensor  Mc 
from  Figure  3.  Agreement  between  data  and  synthetics  indicates  the  success  of  the  moment- 
tensor  inversion.  Residual  time  shifts  are  due  to  3D  propagation  effects  not  removed  by  the 
network  average. 

Transverse-component  data  do  not  constrain  mrr  (since  the  corresponding  Green’s  function, 
grr(t),  is  zero  for  a  spherically  symmetric  earth).  Furthermore,  g66(t)  =  -£^(0>  which  results 
in  mee  =  -m^.  The  LHT  inversion,  therefore,  always  recovers  a  double-couple,  vertical  strike- 
slip  source.  However,  since  mrr  is  unconstrained  it  is  possible  to  fit  the  LHT  data  equally  well 
with  an  alternative  source  mechanism,  the  one  corresponding  to  a  45°  dip-slip  fault  with  a  strike 
rotated  by  45°  and  a  scalar  moment  twice  that  of  the  strike-slip  source  [e.g.,  Aki  and  Richards, 
1980].  To  resolve  this  ambiguity  in  MD(m)  and  to  recover  M;((U)  we  apply  our  inversion 
procedure  to  the  body-wave  portion  between  P  and  the  minor-arc  Rayleigh  waves  (R{)  on 
vertical-component  recordings.  We  calculate  synthetic  test  data  using  the  S12_WM13  3D  earth 
model  and  a  moment  tensor  M  =  Mj  I  +  Mq  M/>,  with  Mq/Mj  =  1  and  from  above.  The 
inversion  allows  for  a  non-zero  isotropic  component;  however,  mr6  and  mr(p  are  projected  out. 
Figure  5  shows  M(co)  between  10  and  45  mHz.  Displaying  the  source  mechanism  unit  vector, 
i,  together  with  the  canonical  double  couple,  CLVD,  and  isotropic  vectors  makes  it  particularly 
easy  to  assess  the  contribution  of  the  purely  dilatational  component.  The  angular  distance 
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between  X  and  i  measures  the  deviation  from  an  isotropic  source.  This  is  the  first  time  that  the 
Riedesel  and  Jordan  [1989]  formalism  has  been  applied  to  assess  and  display  significantly  non- 
deviatoric  mechanisms.  The  agreement  between  different  frequency  bands  and  with  the  known 
solution  (lower  right),  although  not  as  good  as  for  the  transverse-component  data,  is  still  rather 
remarkable.  These  synthetic  tests  indicate  that  we  can  recover  both  deviatoric  and  isotropic 
source  mechanisms  from  vertical-component  and  transverse-component  data. 


A 


30  -  35  mHz 


Figure  5.  Mechanism  plots 
(A)  and  P- wave  radiation 
patterns  (B)  from  the  (phase- 
aligned)  moment-tensor  in¬ 
version  using  the  body-wave 
portion  of  vertical-component 
records  arriving  prior  to  R\. 
Uncertainties  in  the  mecha¬ 
nism  vectors  are  depicted  by 
95%  confidence  ellipses 
(heavy  lines).  The  projected 
input  moment  tensor  (ground 
truth)  consists  of  both 
isotropic,  MfL,  and  devia¬ 
toric,  MdMd,  components 
( Md/Mj  =  1)  and  is  dis¬ 
played  on  the  bottom  right. 
Data  and  synthetics  are  calcu¬ 
lated  for  the  same  models  as 
in  Figure  1.  The  inversion  is 
performed  with  mrQ  and  mr(p 
components  projected  out; 
the  moment  tensor  is  no 
longer  constrained  to  be 
deviatoric. 


Lop  Nor  Nuclear  Test  from  92/5/21.  Having  established  the  performance  of  our  moment-tensor 
recovery  methodology  on  a  synthetic  example,  we  next  apply  our  algorithms  to  data  from  the 
Chinese  nuclear  test  conducted  on  92/5/21.  While  this  is  the  largest  Chinese  nuclear  test  to  date, 
the  station  coverage  provided  by  the  IRIS  and  Geoscope  networks,  is  sparse  by  today's  standards. 
Love  waves  are  clearly  identifiable  on  the  LHT  recordings  from  this  event  and  are  of  equal  or 
larger  amplitude  than  the  Rayleigh  waves  (Figure  6),  demonstrating  that  a  significant  component 
of  tectonic  release  is  present. 
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Figure  6.  Vertical-component  Rayleigh-wave  and  transverse  component  Love-wave 
seismograms  for  six  stations  at  varying  epicentral  distance.  A,  and  source-receiver  azimuth,  <j>, 
from  the  92/5/21  Chinese  nuclear  test  near  Lop  Nor.  Seismograms  are  band-pass  filtered  with  a 
six-pole  Butterworth  filter  (comers  at  8.5  and  31.5  mHz).  Traces  are  aligned  on  the  theoretical 
Love  wave  arrival  time;  maximum  trace  amplitude  is  displayed  below  the  individual  traces. 


11 


Using  the  same  analysis  procedure  as  above,  we  invert  these  data  traces  for  the  deviatoric 
moment  tensor.  However,  we  utilize  our  knowledge  of  long-wavelength  aspherical  earth 
structure  to  calculate  3-D  Green's  functions  using  PREM  together  with  the  3-D  model 
S12_WM13.  Figure  7  displays  the  inversion  results  using  transverse-component,  first-orbit  Love 
waves  (Gj).  The  resulting  moment  tensor,  Mplto),  is  a  strike-slip  mechanism  which  shows 
little  variation  between  frequency  bands,  with  the  strike  of  the  faults  plane  varying  only  by  a  few 
degrees. 


Figure  7.  Mechanism  plots  (A)  and  P- wave  radiation  patterns  (B)  from  the  (phase-aligned) 
Love- wave  inversion  for  the  92/5/21  nuclear  test  near  Lop  Nor,  China.  The  analysis  procedure  is 
the  same  as  for  the  synthetic  test  case  described  earlier,  except  that  we  calculate  Green's  function 
seismograms  for  the  aspherical  model  S12_WM13.  The  resulting  moment  tensor,  M0(a)) 
remains  unchanged  across  the  entire  frequency  band. 


This  event  also  generated  Sh  waves  clearly  observable  at  several  stations  at  regional  and 
teleseismic  distances.  Inverting  these  shear  waves  for  Mo(0))  yields  results  consistent  with  the 
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Love-wave  inversion  (Figure  8).  Since  the  S  waves  are  less  influenced  by  strong,  near-surface 


Figure  8.  Mechanism  plots  (A)  and  P-wave  radiation  patterns  (B)  from  the  inversion  of  phase- 
aligned  Sh  waveforms  for  the  92/5/21  nuclear  explosion  near  Lop  Nor,  China.  The  resulting 
moment  tensor,  M0(<u),  is  nearly  identical  to  that  obtained  from  the  Love- wave  inversion. 


Figure  9  depicts  the  source  mechanism  averaged  over  all  frequency  bands  from  the  Gj  and 
Also  shown  is  the  dip-slip  mechanism  that  satisfies  the  data  equally  well.  The 
moments  of  these  two  mechanisms  are,  =  0.85 ±0.05 x  1017  Nm  and 
2MSq  =  1.7±0.1xl017  Nm,  respectively. 
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Figure  9.  Mechanism  plots  (A)  and  P- wave 
radiation  patterns  (B)  averaged  over  all 
frequency  bands  from  the  Love-wave  and  Sh 
inversions.  Due  to  the  ambiguity  intrinsic  in 
inverting  transverse-component  data,  both 
vertical  strike-slip  and  dip-slip  mechanisms 
satisfy  the  data  equally  well. 


The  network-averaged  amplitude  and  time-shift  spectra  for  the  Sh  waves  using  the  average 
moment  tensor  are  shown  in  Figure  10.  While  the  absolute  amplitude  is  uncertain  within  a  factor 
two,  owing  to  the  ambiguity  between  strike-slip  and  dip-slip  mechanisms,  the  amplitude  spectral 
roll-off  suggests  a  decrease  in  the  size  of  the  tectonic  component,  as  long  as  the  source 
mechanism  does  not  change  from  dip-slip  to  strike-slip  between  0.01  and  0.08  Hz.  The  time 
shifts  at  low  frequencies  (At(co)  =  4-7  s)  indicate  a  centroid  time  more  appropriate  for  an 
earthquake  than  an  explosion. 


Figure  10.  Network-averaged  amplitude  ratios,  A(cq)/a(cq),  and  time-shifts,  A t(co),  (triangles 
with  1-c  standard  errors)  measured  from  narrow- band  cross-correlations  of  Sh  waveforms  using 
an  average  moment  tensor  from  Figure  9.  The  solid  lines  indicate  the  data  fit  obtained  from 
the  source- time  function  displayed  in  Figure  1 1. 


Figure  11  depicts  the  source-time  function  associated  with  the  deviatoric  mechanism,  fo(t), 
obtained  from  a  constrained  inversion  of  the  network-averaged  amplitude  and  time-shift  spectra. 
The  inversion  is  performed  using  a  quadratic  programming  approach  that  allows  for  smoothness 
constraints  and  one-sided  constraints  [Ihmle  et  al.,  1993].  To  limit  the  size  of  the  solution  space, 
we  constrain  //>(?)  to  be  non-negative.  The  source-time  function  shows  some  complexity,  with 
a  sharp  pulse  in  moment  release  followed  by  a  long,  smooth  tail  (Figure  2f).  Although  we  must 
still  examine  potential  sources  of  bias,  we  hypothesize  that  the  pulse-like  component  is  caused  by 
scattering  of  energy  from  the  explosion  into  S//-polarized  waves,  perhaps  by  near-source  lateral 
structure,  and  the  smoothly  varying  component  is  earthquake-like  radiation  due  to  tectonic 
release.  If  correct,  this  interpretation  would  favor  a  tectonic  release  model  where  the  explosion 
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triggers  slip  at  nearby  faults  [Aki  and  Tsai,  1972]  over  the  shatter-zone  model  advocated  by 
Archambeau  [1972]. 

Figure  11.  Source-time  function,  fD{t), 
obtained  from  a  constrained  inversion  of  the 
network-averaged  amplitudes  and  time- 
shifts  depicted  in  Figure  10.  Although  not 
unique,  the  source-time  function  provides  a 
good  spectral  fit.  The  pulse-like  component 
is  hypothesized  to  be  due  to  scattering  of 
explosion  energy,  while  the  smoothly 
varying  component,  comprising  about  60% 
of  the  total  moment  can  be  explained  by 
tectonic  release.  The  moment-rate  is  scaled 
assuming  a  dip-slip  mechanism. 
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To  recover  the  complete  moment  tensor  M(t»)  =  M j  (co)  +  M p((o)  and  to  discriminate 
between  the  two  candidate  mechanisms  obtained  from  inverting  the  SH- type  data,  we  use  the 
vertical-component  data.  We  perform  separate  inversion  using  (1)  first-orbit  Rayleigh  waves  and 
(2)  the  body-wave  portion  between  P  and  Rv  The  Rayleigh-wave  inversions  reveal  the  presence 
of  a  strong  isotropic  component  to  the  source  mechanism,  X,  evident  by  the  small  angular 
distance  between  X  and  i  (Figure  12).  With  a  few  exceptions  the  inversions  are  stable  across 
the  entire  frequency  range.  The  mechanism  uncertainties,  quantified  by  the  95%  confidence 
ellipses,  illustrate  a  trade-off  between  isotropic  and  CLVD  mechanisms.  Rayleigh  waves  are 

A 

unable  to  differentiate  between  an  isotropic  component  and  a  vertical  CLVD,  1',  corresponding 
to  horizontal  extension  and  vertical  compression. 

The  body-wave  data  includes  various  phases  (P,  S,  pP,  pS,  sS,  sP,  etc.),  leaving  the  source  at 
different  take-off  angles,  thus  providing  good  constraints  on  the  source  mechanism.  We  exclude 
higher-mode  surface  waves  that  show  a  strong  contamination  due  to  near-surface  heterogeneities. 
Body-wave  phases  turning  in  the  mantle  are  less  affected  by  shallow  structure  than  surface 
waves,  and  should  therefore  allow  us  to  improve  on  the  Rayleigh-wave  moment-tensor  estimates. 
Indeed,  the  body-wave  inversions  show  consistency  between  M(ct))  determined  in  different 
frequency  bands  and  compare  well  with  the  Rayleigh- wave  results  (Figure  13). 


15 


B 


A  B  A  B 

7  -  9  mHz  15  _  17  mu7 


A 


23  -  25  mHz 


Figure  12.  Mechanism  plots  (A)  and  P-wave  radiation  patterns  (B)  from  the  (phase-aligned) 
vertical-component  Rayleigh- wave  (/?j)  inversion  for  the  92/5/21  nuclear  explosion  near  Lop 
Nor,  China.  The  inversion  attempts  to  recover  M(<y)  =  M/(<y)  +  MpCfo);  however,  mrQ  and 
mr(p  are  projected  out.  The  resulting  moment  tensor,  M(<a),  is  relatively  stable  across  the 
frequency  range  and  exhibits  a  significant  isotropic  component,  evident  by  the  deviation  of  the 
mechanism  vector,  X  from  the  great-circle  of  deviatoric  sources.  The  angular  distance  between 
X  and  i  measures  the  size  of  the  isotropic  component. 


Figure  14  displays  a  frequency-averaged  moment  tensor  using  both  surface-wave  and  body- 
wave  data.  Removing  the  isotropic  component  (Mj  =  3.1±0.3xl017  Nm)  leads  to  the 
deviatoric  moment  tensor  (Mq  =  1.9  ±  0.2  x  1017  Nm)  shown.  The  mechanism  plot  reveals  an 
insignificant  contribution  of  a  CLVD  component;  it  is  virtually  indistinguishable  from  a  double¬ 
couple.  The  best-fitting  double-couple  mechanism  is  displayed  separately  to  facilitate  a 
comparison  with  the  thrust-fault  mechanism  obtained  from  the  inversion  of  the  transverse- 
component  data.  These  two  independently  derived  mechanisms  (and  the  scalar  moments)  are 
quite  similar  (Table  1),  reassuring  us  that  the  inversion  of  the  vertical-component  data 
successfully  separated  isotropic  and  deviatoric  components.  The  tectonic  release  component 
determined  from  our  inversions  is  also  consistent  with  the  limited  amount  of  geologic 
information  available  for  this  area  [Matzko,  1994].  Gao  and  Richards  [1994],  investigating 
earthquakes  near  the  Lop  Nor  nuclear  test  site,  find  a  predominance  of  reverse-faulting  events 
striking  NW  to  SE,  in  good  agreement  with  our  results. 
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Figure  13.  Mechanism  plots  (A)  and 
P- wave  radiation  patterns  (B)  from  the 
(phase-aligned)  moment-tensor  in¬ 
version  using  the  body- wave  portion  of 
vertical-component  records  that  arrive 
prior  to  R\  for  the  92/5/21  nuclear  test 
near  Lop  Nor,  China.  The  resulting 
moment  tensor,  M(<y)  (mrQ  and  mr(p 
are  projected  out)  agrees  well  with  the 
predominant  mechanism  found  from 
the  Rayleigh-wave  inversion. 


Figure  14.  Mechanism  plots  (A)  and 
P- wave  radiation  patterns  (B)  averaged 
over  frequency  bands  from  the  surface - 
wave  and  body-wave  inversions.  The 
proximity  of  the  mechanism  vector  X 
to  i  on  the  focal  sphere  illustrates  the 
strong  isotropic  component.  The 
deviatoric  component  of  the  moment 
tensor  is  nearly  a  double-couple,  and 
differs  only  slightly  from  the  best 
fitting  pure  double-couple  mechanism. 
For  comparison,  the  dip-slip  double¬ 
couple  mechanism  obtained  from 
inverting  the  transverse-component 
data  is  also  depicted.  Double-couple 
source  parameters  are  listed  in  Table  1. 


Table  1.  Source  parameters  of  92/5/21  Lop  Nor  nuclear  test 


Data 

/ — \ 

o 

t— ‘ 

I 

I 

r- 

»-H 

o 

2? 

Md/Mj 

<5t 

At 

05+ 

i?j  and  body  waves 

3.1  ±0.3 

1.9  ±0.2 

0.61  ±0.01 

48° 

98° 

320° 

Gi  and  Sh  waves 

0 

1.7  ±0.1 

45° 

90° 

329° 

+Dip,  8,  slip,  A,  and  strike,  <j)s,  of  best-fitting  double-couple  mechanism. 
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Summary  and  Conclusions 

We  have  developed  a  methodology  for  studying  explosion  sources  that  is  based  on  the 
recovery  of  the  broad-band  moment-rate  tensor,  M(o>),  utilizing  both  body  waves  and  surface 
waves.  We  parameterize  the  source  as  being  composed  of  an  isotropic  and  a  deviatoric 
component,  i.e.,  M(co)  =  M/(<y)  +  This  comprehensive  approach  encompasses  many 

source  parameter  diagnostics  that  have  been  traditionally  used  to  discriminate  nuclear  explosions 
from  chemical  explosions  and  earthquakes  and  has  the  potential  to  provide  new  discrimination 
tools.  Our  formalism  uses  synthetic  seismograms  to  improve  the  localization  of  signal 
measurements  in  both  time  and  frequency  domains.  We  employ  body  waves  and  surface  waves 
from  vertical  component  recordings  to  invert  for  M(<y).  Transverse-component  Sh  and  Love 
waves  provide  independent  estimates  for  MD(<o),  albeit  with  an  intrinsic  ambiguity  between 
strike-slip  and  dip-slip  mechanisms. 

In  our  initial  year  of  funding  we  have  accomplished  the  following  tasks: 

(1)  We  have  developed  procedures  for  inverting  3-component  seismic  data  for  moment-tensor 

representations  of  nuclear  explosions.  These  algorithms  have  the  following  features: 

(a)  they  allow  for  the  facts  that  explosion  sources  are  located  at  very  shallow  depths  by 
projecting  out  the  mrd  and  mr(p  components, 

(b)  they  allow  for  a  source  mechanism  composed  of  both  isotropic  and  deviatoric 
components,  and 

(c)  they  incorporate  3-D  propagation  models. 

(2)  We  have  validated  our  moment-tensor  inversion  algorithms  using  a  synthetic  test  case  with 

a  known  moment-tensor  source,  containing  both  isotropic  and  deviatoric  components. 
The  source  inversions  yield  stable  and  self-consistent  results  across  a  broad  frequency 
range  using  both  body  waves  and  surface  waves  (Figures  1-5). 

(3)  We  have  successfully  applied  graphical  display  tools,  developed  for  assessing  deviatoric 

earthquakes,  to  sources  with  a  significant  dilatational  component.  These  plots  provide  an 
immediate  assessment  of  various  features  of  the  mechanism;  for  example,  the  relative 
importance  of  isotropic  and  deviatoric  (tectonic)  contributions  to  the  moment  release 
(Figures  5, 12-14). 

(4)  We  have  successfully  applied  algorithms  to  invert  source  spectra  for  source  time  functions. 

These  inversions  employ  a  quadratic  programming  algorithm  that  allows  various  types  of 
time-dependent  constraints  to  be  imposed,  including  smoothness  constraints  and  one¬ 
sided  constraints. 
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(5)  To  improve  our  understanding  of  tectonic  release  in  a  variety  of  geologic  settings,  we  have 
focused  on  the  Chinese  nuclear  test  site  near  Lop  Nor,  which  has  received  less  attention 
than  the  U.S.  and  Soviet  test  sites.  Using  recordings  from  the  92/5/21  Chinese  nuclear 
test  we  have  obtained  moment-tensor  solutions  for  both  isotropic  and  deviatoric 
mechanisms.  The  deviatoric  component  of  moment  release  for  this  nuclear  test  is 
significant  (  Mt/Mj  =  0.61  ±  0.01).and  has  a  negligible  non-double-couple  component 
(Figure  14).  The  best-fitting  double-couple  shows  a  reverse-faulting  mechanism  (Table 

I) ,  consistent  with  focal  mechanisms  from  earthquakes  near  the  Lop  Nor  nuclear  test  site. 
Transverse-component  data  (Sh  and  Love  waves),  although  ambiguous  by  themselves, 
provide  an  independent  confirmation  for  MD  derived  from  vertical-component 
recordings  (Figure  14).  The  network-averaged  amplitude  spectrum  rolls  off  above  0.04 
Hz,  indicating  a  decrease  in  the  size  of  the  tectonic  component  (Figure  9)  with  frequency. 
The  source-time  function  associated  with  the  deviatoric  mechanism,  fD(t),  shows  some 
complexity,  with  a  sharp  pulse  in  moment  release  followed  by  a  long,  smooth  tail  (Figure 

II) .  Although  we  must  still  examine  potential  sources  of  bias,  we  hypothesize  that  the 
pulse-like  component  is  caused  by  scattering  of  energy  from  the  explosion  into  Sh- 
polarized  waves,  perhaps  by  near-source  lateral  structure,  and  the  smoothly  varying 
component  is  earthquake-like  radiation  due  to  tectonic  release 

Tectonic  release  can  lead  to  biased  estimates  of  M/(a>).  Furthermore,  the  detection 
threshold  for  a  given  explosion  yield  also  depends  on  the  amount  of  tectonic  energy  released  by 
the  explosion.  It  is  therefore  important  to  better  understand  how  geology  and  explosion 
size/depth,  influence  the  amount  of  tectonic  release  -  in  particular  as  a  function  of  frequency. 
Gaining  new  insights  into  how  tectonic  release  scales  with  frequency  and  explosion  yield/depth 
should  also  lead  to  a  better  understanding  on  how  high-frequency  discriminants  such  as  P/Lg  are 
affected  by  tectonic  release.  At  present,  we  are  addressing  these  questions  by  applying  our 
source  recovery  algorithms  to  data  from  a  number  of  nuclear  tests  of  varying  sizes.  To  better 
quantify  the  importance  of  tectonic  release  as  a  function  of  frequency,  we  are  adapting  our  source 
inversion  tools  to  higher  frequencies.  We  are  also  investigating  formulations  for  directly 
inverting  for  the  relative  contributions  of  tectonic  and  isotropic  components  as  a  function  of 
frequency.  At  higher  frequencies,  apparent  attenuation  due  to  scattering  as  well  as  the  effects  of 
frequency-dependent  intrinsic  attenuation  are  expected  to  dominate  the  source  amplitude 
spectrum,  making  an  inversion  algorithm  based  on  the  relative  magnitude  of  tectonic  and 
explosion  components  desirable. 
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